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'• ABSTRACT 

^ ■ Young planets interact with their parent gas disks through tidal torques. An 

imbalance between inner and outer torques causes bodies of mass ^0.1 Earth 
^ ■ masses (M®) to lose angular momentum and migrate inward rapidly relative to 

^ . the disk; this is known as "Type I" migration. However, protoplanets that grow 

Tij- . to gas giant mass, O(1O^)M0 , open a gap in the disk and are subsequently 

^ . constrained to migrate more slowly, locked into the disk's viscous evolution in 

lO . what is called "Type II" migration. In a young planetary system, both Type I 

^ . and Type II bodies likely coexist; if so, differential migration ought to result in 

p . close encounters when the former originate on orbits exterior to the latter. We 

Q . investigate the resulting dynamics, using two different numerical approaches: an 

il3 . N-body code with dissipative forces added to simulate the effect of the gas disk, 

. and a hybrid code which combines an N-body component with a 1-dimensional 

j> . viscous disk model, treating planet-disk interactions in a more self-consistent 

^ , manner. In both cases, we find that sub-gap-opening bodies have a high like- 

lihood of being resonantly captured when they encounter a gap-opening body. 
A giant planet thus tends to act as a barrier in a protoplanetary disk, collect- 
ing smaller protoplanets outside of its orbit. Such behavior has two important 
implications for giant planet formation: First, for captured protoplanets it miti- 
gates the problem of the migration timescale becoming shorter than the growth 
timescale when Mproto ^ 1^^©- Secondly, it suggests one path to forming systems 
with multiple giant planets: Once the first has formed, it traps (or indeed ac- 
cretes) the future solid core of the second in an exterior mean-motion resonance, 
and so on. The most critical step in giant planet formation may thus be the 
formation of the very first one. 



Subject headings: solar system: formation — planets and satellites: formation 
— planetary systems: protoplanetary disks 
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1. Introduction 

Observational evidence strongly suggests that the formation of giant planets is a ubiq- 
uitous process. To date, well over a hundred extrasolar planets have been discovered by 
radial velocity surveys; in at least eleven cases, multiple planets orbit a single star^. At least 
three of the multi-planet systems have pairs of planets locked in mean-motion resonances 
(MMRs), with a 1:2 (GJ 876b,c; HD 82943b,c) or 1:3 (55 Cnc b,c) period ratio. 

Migration due to planet-disk interaction is the leading candidate for the mechanism 
which produces such dynamical configurations, as well as producing the large fraction of 
"hot Jupiters" , exoplanets on orbits significantly smaller than their likely formation distance 
from their parent star. Indeed, the latter effect was first predicted, and investigated in 
detail, well before it was observed (Goldreich & Tremaine 1980; Lin & Papaloizou 1986). 
More recent work has expanded the focus from planet-disk interaction to "planet-planet- 
disk" dynamics, that is, the behavior of multiple interacting planets embedded in a disk. 
Again, theoretical work anticipated observation: Using hydrodynamic simulations, Kley 
(2000) and Bryden et al. (2000) demonstrated convergent planet migration, thus raising the 
possibility of resonant capture, just prior to the discovery of the first extrasolar MMR, the 
1:2 resonant pair of GJ876b and c. Subsequent work (Masset & Snellgrove 2001; Lee & Peale 
2002; Murray, Paskowitz, & Holman 2002; Nelson & Papaloizou 2002; Kley, Peitz, & Bryden 
2004) looked in more detail on the dynamics of capture and evolution of gas giant planets 
in low-order MMRs. 

The above work focused on the interaction of pairs of gap-opening planets, that is, 
planets which have grown large enough to lock themselves into the viscous evolution of their 
parent disk. For typical protoplanetary disk parameters, this imphes planets of gas giant 
(Jovian) size, O(IO^) Earth masses (M®), rather than ice giant (Neptunian) size, 0(10-*^) 
planets. For initially well-separated planets, convergent migration then requires the removal 
of the intervening gas, thus bringing the planets to share a common gap. It is possible, 
however, for giant planet MMRs to be established earlier. Prior to opening a gap, disk tidal 
torques can cause the orbital decay of an embedded body on a timescale much shorter than 
the disk's viscous time. This mode of migration is commonly referred to as Type I, while 
a gap-opening planet co-evolving with its disk is said to undergo Type II migration Ward 
(1997). Hahn & Ward (1996) suggested that the potentially large differences between Type 
I and Type II timescales could lead to bodies of the former class being delivered to, and then 
resonantly trapped by, bodies of the latter class. The concentration of a significant mass of 



^lAU Working Group on Extrasolar planets, http://www.ciw.edu/boss/IAU/div3/wgesp/planets.shtnil; 
California and Carnegie Planet Search, http://www.exoplanets.org 
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protoplanets in exterior MMRs of a pre-existing giant planet could, they argued, promote 
the formation of subsequent giant planets. In this work, we undertake the first part of a 
numerical investigation of this scenario. We make use of two different approaches. In the 
first and simpler one, we take a standard N-body code and add to it forces to mimic the 
effects of planet-disk interaction, guided by analytic work as well as by numerical simulations. 
The second approach is to take this same N-body code and add to it a simple disk model, 
which viscously co-evolves with embedded planets and interacts with them. The disk model 
is one-dimensional; that is, it describes an azimuthally-averaged thin disk, with disk-planet 
torques obtained from analytic formulae rather than from a full calculation of the launching 
of density waves by the planet. This approach is intermediate between the above "N-body 
with disk forces" and an actual 2- (or 3-) dimensional hydrodynamic simulation. With both 
approaches, wc find that for typical protoplanetary disk parameters, single Type I bodies have 
a high likelihood of being resonantly captured by Type II bodies. Wc further demonstrate 
that significant numbers of protoplancts can occupy adjacent resonances, or indeed share a 
single resonance, of a giant planet. This paper is organized as follows: In §2, we give a brief 
summary of the theory of planet-disk interactions. In §3.1, we present simulations performed 
using the simpler of our two approaches, while in §3.2 we repeat these simulations using our 
hybrid "N-body-|- viscous disk" code. We discuss our results in §4, and summarize them in 
§5. 



2. Planet-disk torques, migration and eccentricity evolution 

Here, we provide a brief summary of the theory of planet-disk interaction, focusing on 
the results we will make use of in this paper. A body orbiting in a gas disk launches density 
waves at Lindblad resonances, exchanging energy and angular momentum with the disk; the 
disk thus modifies the body's orbit. The net effect on the body can be obtained by summing 
over all resonances (Goldreich & Tremaine 1980; see also e.g. Goldreich & Sari 2003). One 
begins by expanding the gravitational potential of the perturbing body in a double Fourier 
series. The /, m component has a pattern speed 

^i,m ^ + (l - 'm)Kp/m (1) 

where Qp and Hp arc the azimuthal and cpicyclic frequencies, respectively, of the planet. 
There is a Lindblad resonance between this potential component and disk material at a 
radius r in the disk where 

^{r) ± ^ = n^^^. (2) 
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Also, a corotation resonance occurs where 



l,m 



(3) 



The radial pressure gradient in the disk causes both the epicychc and the azimuthal frequency 
to differ shghtly from their Keplerian values: 



rp dr 



and 



(4) 



(5) 



In the case of a planet on a circular orbit, all components with I ^ m are zero; ^m,m — 
Qp, so all components have the pattern speed of the planet. Furthermore, there is no torque 
contribution at all from corotation resonances. We will focus on this simpler situation for the 
moment. In summing the effect of all Lindblad resonances, one can define a torque density, 
or torque per unit disk radius. Following Ward (1997), the torque density experienced by 
the disk due to a planet of mass M = /iM* on a circular orbit about a primary of mass 
at radius rp, is 



^1 
dr 
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where { = mCs/rK, Cg is the gas sound speed, and ^ is the dimcnsionless satellite forcing 
function. 
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with P = r/vp and 
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(8) 



yi-2/3cos^ + /?2 

being the Laplace coefficient. Following Menou & Goodman (2004), we account in a simple, 
approximate way for the finite thickness of the disk by replacing the Laplace coefficient with 
a softened approximation. 



(9) 




with Kq the modified Bessel function of the second kind, order 0, H ^ Cg/Q the disk scale 
height, and B a dimcnsionless scaling factor. In going from a summing of torques at discrete 
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resonances to a torque density, the wavenumber m is turned into a continuous function of 
radius: 



m(r) 



1/2 



(10) 



Ward (1988), Artymowicz (1993) and Papaloizou & Larwood (2000) study planet-disk 
torques from an eccentric planet. In particular, Papaloizou & Larwood (2000) look at the 
case where eccentricity exceeds disk aspect ratio H/r. They sum over enough / ^ m torques 
for convergence, in order to obtain the net effect on the planet's orbit. However, they neglect 
corotation torques, as these are zero for their model surface density profile (S oc r~^/^). They 
obtain approximate fits for the timescales of angular momentum and eccentricity change: 
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(12) 



In the limit of low planetary eccentricities, these results are broadly consistent with the 
earher work of Ward (1997) and Artymowicz (1993). Most noteworthy is the short orbital 
decay timescale; for constant eccentricity, a/a — tm/'^, so from Eq. 11, it is of order 10^ years 
for a body of mass 10 in a typical model disk. This would seem to present a significant 
problem for the core-accretion model of giant planet formation, since the assembly of the 
solid core likely takes a million years or more (e.g. Lissauer 1987, Kokubo & Ida 2000, 
Thommes, Duncan, & Lcvison 2003, Inaba, Wetherill, & Ikoma 2003), with perhaps several 
million more required to reach runaway gas accretion (itself a very fast process) and acquire 
a Jovian atmosphere (Pollack et al. 1996). 

If a planet does manage to grow to something approaching gas giant mass before plung- 
ing into its parent star, the problem of rapid inward migration becomes less severe. Torque 
increases with planet mass; for a sufficiently massive body, ~ lO^M® for typical disk models, 
the torques overcome disk viscosity and an annular gap about the orbit is opened in the gas 
(Goldreich & Tremaine 1980), provided that the Hill (Roche) radius also exceeds the disk 
scale height (Lin & Papaloizou 1986, 1993). A planet opening a clean gap prevents disk 
material from crossing its orbit, and thus becomes locked into the disk's viscous evolution. 
In general, this entails a sharp decrease in the rate of orbital decay. 



The eccentricity evolution of a gap-opening planet is as yet not well understood. In 
the gapless case, the dominant effect on the planet's eccentricity is damping by coorbital 
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Lindblad resonances. When a gap forms and empties or strongly depletes the coorbital re- 
gion, the evolution of the eccentricity comes down to a delicate balance between the effect 
of non-coorbital eccentric Lindblad resonances (driving) and corotation resonances (damp- 
ing). Though the latter are slightly stronger than the former, Goldreich & Sari (2003) and 
Ogilvie & Lubow (2003) argue that corotation resonances in fact saturate, thus favouring 
net eccentricity excitation. However, Goldreich & Sari (2003) suggest that damping may be 
reestablished once a planet's eccentricity becomes so high that its radial excursion exceeds 
the gap width. In contrast, numerical simulations tend to show only damping of eccentric- 
ities (e.g. Kley, Peitz, & Bryden 2004) unless the planet mass approaches that of a brown 
dwarf (Papaloizou, Nelson, & Masset 2001). 



For the purpose of the subsequent numerical experiments we will assume that, the 
above difficulties notwithstanding, a protostellar disk does manage to form one initial gas 
giant planet, massive enough to open a gap and orbitally co-evolve with its parent disk. 
Taking this as our initial condition, we then study the evolution of sub-gap-opening bodies 
originating at larger radii than the gas giant. To better assess the robustness of our results, 
we repeat the investigation with two different numerical approaches. 



Simulations are performed with a code based on SyMBA (Symplectic Massive Body 
Algorithm), an N-body program optimized for near-Keplerian systems (Duncan, Levison, 
& Lee 1998). SyMBA is based on the symplectic orbit-integration algorithm of Wisdom & 
Holman (1991), with its high efficiency and good energy-conservation properties, and adds 
an adaptive timestep in order to resolve close encounters between bodies. Using a model gas 
disk, we then calculate the angular momentum and eccentricity decay rate for each sub-gap- 
opening body from Eqs. 11 and 12 each timestep, and apply the corresponding additional 
acceleration a — a^(f) + a^r, where 



3. 



Numerical simulations 



3.1. N-body -|- dissipative disk forces 
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The giant planet in our simulations is treated differently, however. For simplicity, we prevent 
it from migrating altogether. This is done in a way that mimics, rather simplistically, the 
effect of being locked in a disk gap: We choose a radius rgap for the center of the gap and a 
fractional gap "width" w, then apply an azimuthal acceleration a'^ as follows: 



/ Q- ''"gap /I r\ 

\a -rgapl +wa 

Since is negative, a'^ is negative at a > rgap, positive at a < rgap, and falls off smoothly to 
zero as a — > rgap. Therefore, if our "gap-opening" planet is displaced from the gap midpoint, 
it experiences a restoring torque towards a — rgap. 



3.1.1. Single growing protoplanet 

We begin with a simple case: A 300 (approximately Jupiter-mass) body starts on a 
circular orbit at 5 AU about a 1 Mq star, centered in a "disk gap" as described above. We 
place a second body of mass 1 on a circular orbit at 15 AU. A model disk with surface 
density 

Sg = 1000 (^)"'gcm-^ (16) 

and scale height 

/ ^ N 5/4 

if = 0.03(— ) . (17) 

is used to calculate dJ/dt and de/dt. After 5 x 10^ years of simulation time, the outer 
body's mass is increased linearly in time to 30 over the next million years. The resulting 
orbital evolution is shown in Fig. 1. At its original mass of 1 M®, the outer body migrates 
inward until it is captured in the larger body's 1:2 exterior mean motion resonance. Once 
its mass has reached M 2.5Mq, at i 5.5 Myrs, it jumps from the 1:2 to the 2:3 MMR. 
It remains there through the rest of its growth and until the end of the simulation at 5 
Myrs, showing no signs of orbital instability. We perform other simulations with different 
migration rates; we achieve this by simply scaling the disk surface density uniformly, while 
keeping the disk power law index and scale height profile the same (note that J/Jcc if^S~^ 
for small eccentricity). Some of these other results arc also shown in Fig. 1. When E(lAU) = 
1500 gcm~^, the growing body undergoes a second resonance jump, from the 2:3 to the 3:4 
MMR. When S(lAU) = 6000 gcm~^, the body does the same, then also leaves the 3:4 
resonance, but no closer resonance is able to capture it, and it gets all the way to the gap- 
opener's orbit (and in the case pictured, physically collides with it). Our experiment thus 
suggests that, for parameters relevant to a protoplanetary disk — T,g ~ (10^ — 10^)gcm'~^, 
H/r ~ (10~^ — 10~^)r — single bodies of gas giant core/ice giant mass, ~ (10° — 10^) M^, have 
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a high hkehhood of being captured upon encountering a first-order (n:n-|-l) mean-motion 
resonance of a gas giant-sized body. 



3.1.2. Multiple protoplanets 

We now turn to the question of what happens when multiple sub-gap-opening bodies 
encounter a gas giant. Dropping several fully-formed Neptune-mass planets into a disk makes 
for a rather questionable initial condition; as in the previous section, we therefore begin with 
small bodies whose masses we increase over time. This time, we do so in a manner which, 
at least in a broad qualitative sense, conforms to how planet formation proceeds. We make 
use of the prescription for post-runaway or "oligarchic" (Kokubo & Ida 1998) growth in 
the form given by Thommes, Duncan, & Levison (2003). This describes the mass accretion 
rate of a mass Mpp protoplanet growing together with an oligarchy of similar neighbours, all 
embedded in a disk of much smaller planetesimals of characteristic mass nipismi- 

dM„ 



— ^ ^ AM^/' (E, - BM^/^) , (18) 

CLL 

where is the initial local surface density of the planetesimal disk, and A and B are given 

by 

~ 4/15 1/3 1/10 2/15 ' — L 2 vl9j 
PplsmlPpP ^PP "^plsml PP 

Cd is the drag coefficient, a dimensionless constant ~ 1; fe is the spacing in Hill radii between 
adjacent protoplanets; is the mass of the central star; pg k, T,g/ H is the (volume) density 
of gas; and ppp and ppismi are the densities of protoplanet and planetesimals, respectively. 

We put protoplanets, initially of mass 0.1 M®, on circular orbits between 8 and 20 
AU. Their masses are changed over time according to Eq. 18, with multiplied by a 
randomly-generated factor between 0.25 and 1.75 in order to mimic the stochastic nature 
of the accretion. E^ and H are given by Eqs. 16 and 17. We use h = 10, Ppp = Ppismi = 

^plsml 



1.5 gem ^, rupismi = 10 M^; the latter two amount to a planetesimal size of ~ 0.01 km. 



A typical example of the orbital evolution which results is shown in Fig. 2. The outer 
bodies accelerate their inward migration as they grow in mass, and one after another is reso- 
nantly captured by the "Jupiter" at 5 AU. As they continue to grow and the applied torque 
increases, some bodies jump to a closer-in resonance, analogous to the behavior seen in Fig. 
1. At the end of the simulation at 5 Myrs, five bodies share the 2:3 MMR with the gas giant. 
A sixth body is further outward; examination of the different possible resonance angles (e.g. 
Murray & Dermott 1999) shows it to be in a 4:7 MMR with the gas giant, and simultane- 
ously, in a 6:7 MMR with the bodies which are themselves in the 3:2 resonance with the 
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106 2x106 3x106 4x106 5x106 



t (yrs) 

Fig. 1. — Simulation of growing, sub-gap-opening (Type I) body in a gas disk, using the 
"N-body -|- dissipative disk forces" scheme. The simulated disk has a surface density of 
^lAul'^/AU)"^ gcm~^. In the top panel, the orbital evolution is shown for three cases: Siau = 
1000 gcm^^ (black), Siau = 1500 gcm~^ (green), and Siau = 6000 gcm~^ (red). In each 
case the body's semimajor axis a, pericenter q = a(l — e) and apocenter Q = a(l + e) are 
plotted. A body of mass 300 orbits at 5 AU; the Type I body spends time trapped in 
several of its mean- motion resonances (blue dashed lines). The time evolution of the Type 
I body's mass is the same in each of the three cases, and is shown in the bottom panel. 
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giant (4:7 = 2:3x6:7). Meanwhile, all the 2:3 bodies are in a 1:1 commensurabihty with 
each other. This seemingly precarious configuration is attainable here because the tendency 
of the bodies to gravitationally scatter each other as they approach to such close range, is 
kept in check by the strong damping of eccentricities which is simultaneously applied. Other 
outcomes are shown in Fig. 3. The 2:3 resonance is strongly favoured as the final location 
of the growing protoplanets, with the smallest bodies tending to end up further out, in a 5:9 
(5:6) or 4:7 (6:7) MMR with the gas giant (the 2:3 MMR bodies), respectively 

3.2. Hybrid N-body + viscous disk 

We now repeat our experiments with a different numerical approach. The code used here 
is called VDisk-SyMBA; as the name suggests, it also uses SyMBA for evolving planetary 
bodies, but the disk is treated in a more sophisticated way: We co-evolve with the planets 
a one-dimensional disk model, which not only changes in time due to its own viscosity, but 
also exchanges angular momentum with the N-body part, modifying its own density profile 
as well as the orbits of the planets. A similar computational approach, though without the 
N-body component, was utilized by Lin & Papaloizou (1986) to follow the evolution of single 
protoplanets in a disk. Since radius is its only dimension, the disk is of necessity vertically 
and azimuthally averaged. However, both vertical and azimuthal structure are implicitly 
included in the model. For the former, we simply assign a scale height; for the latter, 
planet-disk interaction — a process which after all works by raising azimuthally asymmetric 
perturbations in the disk — is handled by way of the torque density formulation of Eq. 
6 applied to the azimuthally averaged disk surface density. This approach to simulating 
planet-disk interactions is thus intermediate in complexity between our first method and a 
full hydrodynamic simulation. The fidehty of the planet-disk interactions is only as good as 
the torque density used, but because the disk is one- dimensional, this code has a significant 
speed advantage over two- or three-dimensional hydrodynamic codes. In practice we find 
that millions of years can typically be simulated in one day of real time on a desktop- 
class machine, thus simulations spanning the entire lifetime of a protoplanetary gas disk, 
constrained observationally to be ~ 10^ years or less (Haisch, Lada, & Lada 2001) are within 
easy reach. 

The gas disk is divided into radial bins of arbitrary size at the start of the simulation. 
The most stringent upper limit on bin size comes from the necessity of properly resolving 
planet-disk torques; accurately tracking the disk's viscous self-evolution is much more for- 
giving in that respect. We generally find it necessary to scale the bin size to be some fraction 
of the local disk scale height. 
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Fig. 2. — A simulation of multiple growing protoplanets outside a giant planet (fixed at 5 
AU) using the "N-body + dissipative disk forces" scheme. Top panels: Orbital evolution of 
semimajor axis, peri- and apocenters of the bodies. The left panel shows the initial migration, 
while the one on the right shows the full simulation with a a smaller vertical scale. Bottom 
left: Evolution of the protoplanet masses by the "oligarchy" prescription. Bottom right: 
Final state showing resonance locations. 
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Fig. 3. — Snapshots taken at 5 x 10^ years of a set of ten simulations (differing only in 
the initial orbital phases of the bodies) of a population of growing protoplanets on orbits 
exterior to a gap-opening gas giant. The locations of some of the giant's exterior mean- 
motion resonance centers are shown (vertical lines). Most of the surviving protoplanets are 
locked in the 2:3 resonance, typically with one additional body in the 4:7 and 5:9 resonance. 
In one case (G), there is a single much larger protoplanet, formed from the merger of several 
smaller ones, in the 3:4 resonance. 
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The governing equation for the surface density of a Keplerian viscous disk is 

(20) 



dE _1 d 
dt r dr 



' 1/2 -9 , ^ 1/2, rV2 dT 



where dT/dr is the torque density in the disk due to the planet (Eq. 6). For the viscous 
self-evolution of the disk, we solve this equation each timestep, with the second term set to 
zero, using an explicit upwind differencing scheme. The planet-disk torques are computed 
and applied to both planet and disk in a separate step. An important part of this is to 
compute in each bin; although we assume perfectly Keplerian gas velocities for the disk's 
viscous evolution, the torque asymmetry felt by sub-gap-opening bodies is, as mentioned in 
§2, primarily due to the sub-Kcplerian gas orbital speed shifting resonance locations. We 
calculate Vt in each bin at each timestep using the finite- difference version of Eq. 4. We then 
make the approximation that k = Vt. 

The planet-disk torques are obtained using the expression for torque density, Eq. 6. We 
are thus computing the torques for a circular orbit; when e 7^ 0, we use a simplified version 
of Eq. 12 (without the eccentricity-dependent part in square brackets) to separately change 
the eccentricity. For gas surface density Eg, we average among the value at the semimajor 
axis, pericentre and apocentre of the orbit. This is a somewhat ad-hoc prescription for the 
eccentricity evolution, but it does reproduce the analytically-obtained results for a non-gap- 
opening body, while for a gap-opening planet, damping of eccentricities does not commence 
until the planet's radial excursion causes it to overrun its gap. Thus we have implemented a 
compromise of sorts between numerical results on the one hand (e.g. Papaloizou, Nelson, & 
Masset 2001; Kley, Peitz, & Bryden 2004) and the analytic work of Goldreich & Sari (2003) 
and Ogilvie & Lubow (2003) on the other. It should also be pointed out that implicit in 
our approach is the assumption that angular momentum is deposited in the disk right at 
the location of a resonance. This is a good assumption as long as the torquing planet has 
rn greater than the scale height H of the disk, in this case waves launched at resonances 
will indeed be nonlinear and shock immediately, thus depositing their angular momentum 
right away (Lin & Papaloizou 1993). However, for bodies having rn < H, a more accurate 
treatment would include modelling the nonlocal deposition of angular momentum in the disk 
(Takeuchi, Miyama, Sz Lin 1996; Rafikov 2002). The effect of our simplification is that we 
underestimate the gap-opening threshold mass for bodies with rn < H. 



3.2.1. Single protoplanet 



We generate an initial gas disk of similar surface density and scale height profile as in 
§3.1.1. We adopt a standard a-parameterization for the viscosity of the disk, u — aCgH, 
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with a — 10~^. A 300 M® body is, again, started on a circular orbit at 5 AU, while a 1 
body is started on a circular orbit at 15 AU. As before, after 5 x 10^ years, its mass is 
linearly increased to 30 M® over the next 10^ years. The orbital evolution of both bodies, 
together with the disk, is shown in Fig. 4 

The most obvious difference in the evolution, with respect to Fig. 1, is that the inner 
planet also migrates. It immediately opens a gap in the disk, initially moving outward slightly 
as the nearby disk material rearranges itself about the freshly-opened gap. Thereafter, the 
planet moves inward along with the accretion flow of the disk, going from 5 AU to 2 AU 
in just under 2 Myrs. The outer planet migrates much faster, at a similar rate as for the 
Eg(lAU) = 1000 gcm^^ case in §3.1.1, and is captured into the 1:2 MMR of the inner planet 
after a bit more than 10^ years. Unlike the cases shown in Fig. 1, it remains in the 1:2 
resonance as it grows. The strength with which the outer body is pushed toward the inner 
body is somewhat reduced by the inward migration of the inner body, as well as by the 
time evolution of the gas disk (S^ decreases as the disk spreads). Also, as can be clearly 
seen in Fig. 4, as the outer body becomes more massive (from --^ 8 x 10^ years onward) it 
begins to open a gap of its own. As this happens, gas begins to be trapped between the 
two bodies; the more effective the trapping, the more this gas acts to push the bodies apart, 
counteracting the asymmetric torque felt by the outer body. Another difference with respect 
to Fig. 1 is that the gap-opening body temporarily acquires a substantial eccentricity, such 
that its radial excursion becomes equal to the gap width, where the onset of strong damping 
then halts eccentricity growth. At the same time, the smaller body's eccentricity decreases 
temporarily (examination of the 2:1 resonance angles shows a coincident change in their 
libration centers). This behavior appears to arise in a fairly narrow range of parameter 
space, wherein the torque on the outer body is, for an extended period of time, just below 
the threshold needed to move it from the 1:2 to 2:3 resonance. It is not apparent in any of 
the other simulations presented below, and we defer a closer examination to future work. 

We perform a second simulation with a more massive gas disk. Eg 3000 (r/1 AU)^^ gem" 
The results are shown in Fig. 5; this time, the outer planet goes from the 1:2 to the 2:3 
resonance with the inner almost as soon as its growth commences. This puts it right at the 
outer edge of the gap. 

3.3. Multiple protoplanets 

Using the same prescription as in §3.1.2, we again grow multiple protoplanets originating 
as 0.1 M0 bodies exterior to the gas giant. A typical simulation is shown in Fig. 6. In 
addition, snapshots at 1.5 x 10^ years of a set of similar simulations (differing only in the 
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5.0x105 1.0x106 1.5x105 

t (yrs) 

Fig. 4. — Orbital evolution of an inner gap-opening (Type II) 300 planet and an outer 
sub-gap-opening (Type I), initially 1 body. For each body, the semimajor axis a (black), 
pericenter q = a(l — e) and apocenter Q = a{a + e) (both green) are plotted. The locations 
of the gap-opening body's 1:2 and 2:3 mean-motion resonances (MMRs) are also shown 
(upper and lower yellow dashed line, respectively). The 1 body grows linearly to 30 
between 5 x 10^ yrs and 1.5 x 10^ yrs; it remains trapped in the 1:2 MMR until the end of 
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5.0x105 1.0x106 1.5x105 

t (yrs) 

Fig. 5. — The simulation shown in Fig. 4, repeated with a higher initial gas disk surface 
density of S 3000(r/l AU)~^ gcm~^. This time, the outer body transitions from the 1:2 
to the 2:3 MMR at 5 x 10^ years. 
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randomly-generated initial orbital phases of the bodies) are shown in Fig 7. All surviving 
protoplanets at this time are locked in either the 2:3 or the 1:2 exterior MMR of the gap- 
opening giant, with the majority residing in the former. As in the simulations of §3.1.2, 
multiple protoplanets — as many as four in this case — end up sharing the 2:3 resonance. 

4. Discussion 

The numerical simulations we have performed here suggest that when rapidly-migrating, 
sub-gap-opening protoplanets approach a gap-opening gas giant, resonant capture is a likely 
outcome. The two different approaches we use in our simulations — first, dissipative forces 
added to an N-body code, and then, a more self-consistent treatment of the disk which in- 
cludes gap formation and viscous accretion — both produce similar results. A particularly 
interesting feature common to both sets of simulation results is the "stacking" of multiple 
multi-Earth-mass bodies in a single mean-motion resonance of the gas giant. Such seemingly 
precarious resonant configurations can be assembled in our simulations because the eccen- 
tricity damping we apply to sub-gap bodies stabilizes them against catastrophic mutual 
scattering with pre-existing resonant bodies as they approach a MMR. 

The ability of a giant planet locked into the viscous transport of the disk to arrest the 
orbital decay of 0(10*^ — 10^) bodies, which, at least in a simple disk model like the one 
used here, move inward much faster than the disk material, has important implications for 
giant planet formation. In the classic core-accretion model (Pollack et al. 1996), most of 
a giant planet's "childhood", several Myrs in length, is spent as a ~ 10^ body waiting 
to accrete a massive gas envelope. However, the work of Ward (1997) and e.g. Papaloizou 
& Larwood (2000), summarized in §2, implies an orbital decay timescale for such bodies 
of only ~ 10"^ years. Thus, having cores held up by a more slowly-migrating gas giant can 
potentially be a big help in enabling them to reach core mass at all. However, even a gap- 
opening planet does not offer a safe haven indefinitely. As can be seen in the simulations 
of §3.2, the gas giant undergoes significant migration in only 2 Myrs, moving in more than 
1 AU/Myr on average. And this, it should be stressed, is with a rather low a viscosity 
parameter of 10~^. With a = 10~^, a gap-opening planet starting at 5 AU has well below 1 
Myr before it plunges into the star, barring some sort of "parking mechanism" at the inner 
edge of the disk (e.g. a central magnetospheric cavity; Lin, Bodenhcimcr, & Richardson 
1996). In other words, for a sufficiently large disk viscosity. Type 11 migration becomes 
just as problematically fast as Type I migration. Thus, in order for a gas giant to be of 
use as a safety net for smaller bodies, it must exist in a disk, or a region of a disk, that 
has fairly low viscosity. Just such a possibihty is offered by the layered accretion or "dead 
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t (yrs) 

Fig. 6. — Simulation of the orbital evolution of a population of growing protoplanets outside a 
gap-opening gas giant, plotted as in Figs. 4 and 5. The protoplanets migrate inward rapidly, 
and one protoplanet after another encounters the exterior mean-motion resonances of the 
giant. After a period of vigorous interactions lasting until ~ 8 x 10^ years, the protoplanets 
all occupy either the 1:2 or the 2:3 MMR. 
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Fig. 7. — Snapshots taken at 1.5 x 10^ years of a set of ten simulations (differing only in 
the initial orbital phases of the bodies) of a population of growing protoplanets on orbits 
exterior to a gap-opening gas giant. The locations of the exterior 2:3 and 1:2 mean-motion 
resonances of the gap-opening planet are shown as vertical lines. All surviving protoplanets 
are locked in one or the other resonance. 
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zone" disk model (Gammie 1996; Matsumura & Pudritz 2003), in which the bulk of the 
disk transports angular momentum only at small radii, < 1 AU, and beyond several AU. In 
between, angular momentum is only transported in the disk's surface layers. This is because 
the assumed source of disk viscosity, magnetorotational instability (Balbus & Hawley 1991), 
only operates in a sufficiently ionized disk; at small radii the star can provide this ionization 
and at larger radii cosmic rays take over, but at intermediate radii, the column depth of 
the disk is too high to allow cosmic rays to penetrate all the way through. The result is 
that from less than 1 AU to several AU (the exact bounds depend on the details of the 
model) the disk has a low viscosity, and thus gas — and gap-embedded planets — will accrete 
slowly through this region, giving planet formation more time to play out. Even then, the 
ultimate survival of whatever planets are produced may require the removal of at least part 
of the outer disk by photocvaporation rather than accretion (Shu, Johnstone, & Hollenbach 
1993; Hollenbach, Yorkc, & Johnstone 2000; Matsuyama, Johnstone, & Murray 2003). It 
has been suggested that a giant planet may promote the accretion of further planets due to 
the enhancement of disk surface density near the edges of a freshly-excavated gap (Lin & 
Papaloizou 1979). The outer gap edge, in particular, may locally concentrate solids because 
the transition from super- to sub-Keplerian gas orbital speed there ought to causes a pile-up 
of planetesimals in a narrow radial range (Bryden et al. 2000). We demonstrate that, in 
agreement with the suggestion of Hahn & Ward (1996), an analogous effect exists when 
solids grow to larger size, where disk tides rather than aerodynamic drag are the primary 
agent of gas-solids interaction. Indeed, both effects — the pile-up of small planetesimals and 
the resonant capture of large protoplanets — arc likely to operate in concert. This brings 
up a shortcoming in our prescription for protoplanet growth: We use an estimate based 
on orderly post- runaway ("oligarchic") growth, which assumes a smooth surface density of 
solids and a uniform spacing of neighbouring protoplanets. Both assumptions will certainly 
not hold up well when a local enhancement of solids exists, and especially when multiple 
bodies are competing to accrete in the same mean-motion resonance. However, our intent 
here has been merely to make the accretion prescription plausible in the broadest sense; a 
more detailed study of the accretion of multiple resonantly-locked, co-orbital protoplanets is 
a subject for future work. 

A shortcoming common to both of our simulation schemes is that wc neglect the effect 
of corotation resonances on sub-gap-opening planets. Since their strength is proportional to 
the gradient of the the quantity (fiSg//?^), this is more of an issue for gap-opening planets, 
which have large gradients at the gap edges. These resonances are pivotal in determining 
the rate — indeed the sign — of the eccentricity evolution of planets in gaps (Goldreich & Sari 
2003; Ogilvie & Lubow 2003). In modeling gap-opening planets, we also make no explicit use 
of corotation resonance torques, but rather use a compromise among the competing models 
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for the net effect of all resonances on the eccentricity evolution: Eccentricity is neither driven 
nor damped while the planet is fully in its gap, but is damped once it becomes eccentric 
enough to run into the gap edges. A more general caveat is that our simulation results 
are contingent on the validity of treating planet-disk torques of multiple nearby bodies in 
an azimuthally-averaged, linearly-superposed manner. It will be very interesting to sec to 
what extent full hydrodynamic simulations reproduce the resonant configurations we have 
generated here, especially the extreme case of multiple Type I bodies sharing a Type II 
body's MMR. 

The simulations performed here lend support to the idea that the prc-cxistcncc of a 
giant planet makes the formation of subsequent planets easier. If so, the formation of the 
first giant planet may be the critical event in a protoplanetary disk. But with no "safety 
net" of its own, how can the first sohd core be assembled fast enough? Menou & Goodman 
(2004) offer one possibility: In a more complex disk model than is usually used for such 
cak:Tilations (including the ones in this work), opacity transitions can create local migration 
bottlenecks which stretch the Type I timescale to ~ 10^ years. The work of Cuzzi & Zahnle 
(2004) is also suggestive: They show that icy planetesimals undergoing orbital decay due to 
aerodynamic gas drag can produce and maintain a massive plume of water vapor once they 
reach the water sublimation radius ("snow line"). If much of this water recondenses just 
outside the snow line, a large local enhancement in the surface density of solids may result. 
At that radius, accretion will be more rapid and produce bigger bodies than elsewhere. 
Furthermore, although we have been working from the perspective of the core accretion 
model, one can also imagine a hybrid scenario wherein the first giant planet forms rapidly 
via the competing disk instability model (Cameron 1978, Boss 1997, Mayer et al. 2002), 
thus avoiding the Type I issue altogether, and subsequently facilitates the formation of its 
successors by core accretion. This initial, innermost planet (or planets) may well end up 
migrating into the star, leaving only the successors to populate the mature system. That 
could circumvent some of the inconsistencies with observations commonly associated with 
the disk instability model, namely higher-than-detected planet masses and (in the context 
of our Solar System) insufficiently enhanced planet metallicities. Finally, regardless of the 
formation mechanism, it is worth remembering that only about a quarter of detected multi- 
planet systems contain resonant pairs. Therefore, an important question to be addressed in 
future work is how readily mean-motion resonances can be broken as a system evolves. 
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5. Summciry 

We have performed numerical simulations which suggest that gap-opening (Type II) 
gas giant-sized bodies are very efficient at resonantly capturing sub-gap-opening (Type I) 
bodies of Earth to Neptune size in their exterior mean-motion resonances. In fact, multiple 
such bodies may accumulate in a single resonance. The local concentration of such a large 
amount of mass for significant lengths of time stands in sharp contrast to the problem which 
arises when considering the core-accretion growth of a single giant planet in isolation, namely 
that the growth timescale exceeds the migration timescale long before the critical mass for 
runaway accretion of a gas envelope (~ 10^ M^) is reached. Hence, the first gas giant to form 
may be the most important, paving the way for its successors. Also, these results suggest 
one way to build mean-motion resonant systems of giant planets like those discovered by 
radial-velocity surveys to be orbiting GJ 876, HD 82943, and 55 Cnc: Rather than giant 
planets forming separately and then undergoing convergent migration, resonance locking 
may happen much earlier, when the second planet is still just a growing core. 

I would like to thank G. Bryden, W. Kley, D. Lin, M.-H. Lee, G. Marcy and F. Masset 
for stimulating and informative discussions, the Kavli Institute for Theoretical Physics for 
its hospitality during the early part of this work, and an anonymous referee for valuable 
comments and suggestions. This work is supported by the Natural Sciences and Engineering 
Research Council of Canada. 
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